Optimal design of spatial distribution networks 



Michael T. Gastner^-^ and M. E. J. Newman^-^ 

^ Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501 
^Department of Physics, University of Michigan, Ann Arbor, MI 48109 
^Center for the Study of Complex Systems, University of Michigan, Ann Arbor, MI 48109 

We consider the problem of constructing public facilities, such as hospitals, airports, or malls, 
in a country with a non-uniform population density, such that the average distance from a per- 
son's home to the nearest facility is minimized. Approximate analytic arguments suggest that the 
optimal distribution of facilities should have a density that increases with population density, but 
does so slower than linearly, as the two-thirds power. This result is confirmed numerically for the 
particular case of the United States with recent population data using two independent methods, 
one a straightforward regression analysis, the other based on density dependent map projections. 
We also consider strategies for linking the facilities to form a spatial network, such as a network of 
flights between airports, so that the combined cost of maintenance of and travel on the network is 
minimized. We show specific examples of such optimal networks for the case of the United States. 



I. INTRODUCTION 

Suppose we are given the population density p(r) of a 
country or province, by which we mean the number of 
people per unit area as a function of geographical posi- 
tion r. And suppose we are charged with choosing the 
sites of p facilities, such as hospitals, post offices, su- 
permarkets, gas stations, or schools, so that the mean 
distance to the nearest facility averaged over the popula- 
tion is minimized. In most countries population density 
is highly non-uniform, in which case a uniform distribu- 
tion of facilities would be a poor choice: it gains us little 
to build a lot of facilities in sparsely populated areas. A 
more sensible choice would be to distribute facilities in 
proportion to population density, so that a region with 
twice as many people has twice as many facilities. But 
this distribution too turns out to be suboptimal, because 
we also gain little by having closely spaced facilities in 
the highly populated areas — when facilities are closely 
spaced the typical person is not much further from their 
second-closest facility than from their closest, so one or 
the other can often be removed with little penalty and 
substantial savings. 

As we will see, the ideal solution to this problem lies 
somewhere between these two extremes, with the density 
of facilities increasing as the two-thirds power of popu- 
lation density, a prediction that we verify using simula- 
tions and visualizations based on cartograms, with ac- 
tual population data for the United States. In addition, 
one is often interested in connections between facilities, 
such as flights between airports Q or transmission lines 
between power stations . In the second half of this pa- 
per we generate networks based on a simple model that 
optimizes network topology with respect to the cost of 
maintaining and traveling across the network. Depend- 
ing on the benefit function chosen, we find structures 
ranging from completely decentralized networks to hub- 
and-spoke networks. 



II. OPTIMAL DISTRIBUTION OF FACILITIES 

We wish to distribute p facilities over a two- 
dimensional area A such that the objective function 

/(ri...rp)= / p{r) min |r-ri|d^r, (1) 
Ja ie{i...p} 

is minimized. Here {ri . . . Vp} is the set of positions of 
the facilities and p{r) is the population density within 
the region A of interest. This objective function is pro- 
portional to the mean distance that a person will have to 
travel to reach their nearest facility. 

Seemingly simple, this so-called p-median problem has 
been shown to be NP-hard , so in practice most studies 
rely either on approximate numerical optimization or ap- 
proximate analytic treatments H. A number of different 
approaches have been used 0, fsl M H Q ; the calculation 
given here is essentially that of Gusein-Zade [lof . 

Our p facilities naturally partition the area A into 
Voronoi cells. (The Voronoi cell Vi for the «th facility 
is defined as the set of points that are closer to than 
to any other facility.) Let us define s(r) to be the area 
of the Voronoi cell to which the point r belongs. In two 
dimensions a person living at point r will on average be a 
distance g[s{r)]-^/^ from the nearest facility, where 5 is a 
geometric factor of order 1, whose exact value depends on 
the shape of the Voronoi cell, but which will in any case 
drop out of the final result. The distance to the nearest 
facility averaged over all members of the population is 
proportional to 

f^g f p(r)[s(r)]i/2d\ (2) 

J A 

where we are making an approximation by neglecting 
variation of the geometric factor g between cells. 

The value of s(r) is constrained by the requirement 
that there be p facilities in total. Noting that s(r) is 
constant and equal to s{ri) within Voronoi cell Vi, we see 



FIG. 1: Facility locations determined by simulated anneal- 
ing and the corresponding Voronoi tessellation for p = 5000 
facilities located in the lower 48 United States, based on pop- 
ulation data from the US Census for the year 2000. 
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FIG. 2: Facility density D from Fig.Qversus population den- 
sity p on a log-log plot. A least-squares linear fit to the data 
gives a slope of 0.663 ± 0.002 (solid line). 



that the integral of [s(r)] ^ over Vi is 



[s(r)]-i dV = [s(r,)]-i / = 1. (3) 



Summing over all Vi, we can then express the constraint 
on the number of facilities in the form 



[s{r)] 



-id^r 



(4) 



Subject to this constraint, optimization of the mean 
distance / above gives 
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-0, 
(5) 

where a is a Lagrange multiplier. Performing the func- 
tional derivatives and rearranging for s(r), we find s(r) = 
[2a/gp(r)]2/3. The La grange multiplier can be evaluated 
by substituting into Eq. and we arrive at the result 



^ ' s(r) ^/[p(r)]2/3d2r' 



(6) 



where we have introduced the notation D{r) = [s(r)]~^ 
for the density of the facilities. 

Thus, if facilities are distributed optimally for the given 
population distribution, their density increases with pop- 
ulation density but does so slower than linearly, namely 
as a power law with exponent | 27]. This density places 
most facilities in the densely populated areas where 
most people live while still providing reasonable service 
to those in sparsely populated areas where a strictly 
population-proportional allocation might leave inhabi- 
tants with little or nothing. 

The derivation above makes two approximations: it 
assumes that the geometric factor g is the same for all 



Voronoi cells and that s(r) is a continuous function. Nei- 
ther assumption is strictly true, but we expect them to 
be approximately valid if p varies little over the typi- 
cal size of a Voronoi cell. As a test of these assump- 
tions, we have optimized numerically the distribution of 
p = 5000 facilities over the lower 48 states of the United 
States (Fig. using population data from the most re- 
cent US Census jllj . which counts the number of resi- 
dents within more than 8 million blocks across the study 
region. To create a continuous density function p, we 
convolved these data with a normalized Gaussian dis- 
tribution of width 20 km. The facility locations were 
then determined by optimizing the full p-median objec- 
tive function by simulated annealing (12i] . 

The relation D oc p^/a ^g^j^ |-,g tested as follows. First, 
we determine the Voronoi cell around each facility. Then 
we calculate D{r) as the inverse of the area of the cor- 
responding cell and p as the number of people living in 
the cell divided by its area. Figure [21 shows a scatter plot 
of the resulting data on doubly- logarithmic scales. If the 
anticipated |-power relation holds, we expect the data to 
fall along a line of slope |. And indeed a least-squares fit 
(solid line in the figure) yields a slope 0.663(2), in good 
agreement with the theoretical prediction. 

Some statistical concerns might be raised about this 
method. First, we used the Voronoi cell area to calculate 
both D and p, so the measurements of x- and y- values in 
the plot are not independent, and one might argue that 
a positive slope could thus be a result of artificial corre- 
lations between the values rather than a real result [l3l |. 
Second, it is known that estimating the exponent of a 
power law such as Eq. (jSJ from a log-log plot can in- 
troduce systematic biases 0|. In the next section, we 
introduce an entirely different test of Eq. ((HJ that, in ad- 
dition to being of interest in its own right, suffers from 
neither of these problems. 
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III. DENSITY-EQUALIZING PROJECTIONS 

If we neglect finite-size effects, it is straiglitforwarcl to 
demonstrate tfiat optimally located facilities in a uni- 
formly populated space lie on the vertices of a regular 
triangular lattice [13 ■ It has been conjectured that for 
a non-uniform population there is a general class of map 
projections that will transform the pattern of facilities 
to a similarly regular structure 163 . The obvious candi- 
date projections are population density equalizing maps 
or cartograms, i.e., maps in which the sizes of geographic 
regions are proportional to the populations of those re- 
gions 

in III in 113 ■ Densely populated regions appear 
larger on a cartogram than on an equal-area map such 
as Fig. n and the opposite is true for sparsely populated 
regions. Since most facilities are located where the pop- 
ulation density is high, a cartogram projection will effec- 
tively reduce the facility density in populated areas and 
increase it where the population density is low. There- 
fore, one might expect that a cartogram leads to a more 
uniform facility density than that shown in Fig. ^ And 
indeed some authors have used population density equal- 
izin g pr ojections as the basis for facility location meth- 
ods |2ll|23. 

In Fig. we show the facilities of Fig. on a pop- 
ulation density equalizing cartogram created using the 
diffusion-based technique of Although the popula- 
tion density is now equal everywhere, the facility density 
is obviously far from uniform. A comparison between 
Fig. n and Ot- reveals that we have overshot the mark 
since the facilities are now concentrated in areas where 
there are few in actual space. 

Equation ((BJ makes clear what is wrong with this ap- 
proach. Since D grows slower than linearly with p, a 
projection that equalizes p will necessarily overcorrect 
the density of facilities. On the other hand, based on 
our earlier result, we would expect a projection equal- 
izing p"^/^ instead of p to spread out the facilities ap- 
proximately uniformly. Hence, one way to determine the 
actual exponent for the density of facilities is to create 
cartograms that equalize p^, x > 0, and find the value 
of X that minimizes the variation of the Voronoi cell sizes 
on the cartogram. This approach does not suffer from 
the shortcomings of our previous method based on the 
doubly-logarithmic plot in Fig. |21 since we neither use 
the Voronoi cells to calculate the population density nor 
take logarithms. One might argue that the Voronoi cells 
on the cartogram are not equal to the projections of the 
Voronoi cells in actual space, which is true — the cells gen- 
erally will not even remain polygons under the cartogram 
transformation. The difference, however, is small if the 
density does not vary much between neighboring facili- 
ties. 

In Fig. 21 we show the measured coefficient of variation 
(i.e., the ratio of the standard deviation to the mean) for 
Voronoi cell sizes on p^ cartograms as a function of the 
exponent x (solid curve). As the figure shows, the mini- 
mum is indeed attained at or close to the predicted value 



of X = |. Fi gurelJh shows the corresponding cartogram 
for this exponent. This projection finds a considerably 
better compromise between regions of high and low pop- 
ulation density than either Fig.^orl^li. 

For comparison, wc have also made the same measure- 
ment for 5000 points distributed randomly in proportion 
to population. Since the density of these points is by 
definition equal to p, we expect the minimum standard 
deviation of the cell areas to occur on a cartogram with 
X = 1. Our numerical results for this case (dotted curve 
in Fig. 01) agree well with this prediction. Comparing 
the solid and the dotted curves in the plot, we see that 
not only the positions of the minima differ, but also the 
minimal values themselves. The lower standard devia- 
tion for the p-median distribution indicates that opti- 
mally located facilities are not randomly distributed with 
a density oc p^^^. Instead, the optimally located facili- 
ties occupy space in a relatively regular fashion reminis- 
cent of the triangular lattice of the uniform population 
case [Hill. 



IV. OPTIMAL NETWORKS OF FACILITIES 

In many cases of practical interest, finding the opti- 
mal location of facilities is only half the problem. Often 
facilities are interconnected forming networks, such as 
airports connected by fiights or warehouses connected by 
truck deliveries. In these cases, one would also like to find 
the best way to connect the facilities so as to optimize 
the performance of the system as a whole. 

Consider then a situation in which our facilities form 
the nodes or vertices of a network and connections be- 
tween them form the edges. The efficiency of this net- 
work, as we will consider it here, depends on two factors. 
On the one hand, the smaller the sum of the lengths of 
all edges, the cheaper the network is to construct and 
maintain. On the other hand, the shorter the distances 
through the network between vertices, the faster the net- 
work can perform its intended function (e.g., transporta- 
tion of passengers between nodes or distribution of mail 
or cargo). These two objectives generally oppose each 
other: a network with few and short connections will 
not provide many direct links between distant points and 
paths through the network will tend to be circuitous, 
while a network with a large number of direct links is 
usually expensive to build and operate. The optimal so- 
lution lies somewhere between these extremes. 

Let us define kj to be the shortest geographic distance 
between two vertices i and j measured along the edges 
in the network. If there is no path between i and j, we 
formally set Uj = 00. Introducing the adjacency matrix 
A with elements Aij = 1 if there is an edge between 
i and j and Aij = otherwise, we can write the total 
length of all edges as 

T = J2A,,k,. (7) 

i<j 





FIG. 3: Near-optimal facility location on (a) a cartogram equalizing the population density p and (b) a cartogram equalizing p 
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FIG. 4: The coefficient of variation (i.e., the ratio of the stan- 
dard deviation to the mean) for Voronoi cell areas as they 
appear on a cartogram, against the exponent x of the under- 
lying density p^ . 



We assume this quantity to be proportional to the cost 
of maintaining the network. Clearly this assumption is 
only approximately correct; networked systems in the 
real world will have many factors affecting their mainte- 
nance costs that are not accounted for here. It is however 
the obvious first assumption to make and, as we will see, 
can provide us with good insight about network struc- 
ture. 

The typical cost of shipping a commodity or traveling 
through the network depends on the distances kj as well 
as the amount of traffic Wij (e.g., weight of cargo, number 
of passengers, etc.) that flows between vertices i and j. 
In a spirit similar to our assumption about maintenance 
costs, we assume that the total travel cost is given by 



Z — ^ Wijkj. 

i<j 



(8) 



populations in the Voronoi cells Vi and Vj around i and j, 
so that 



Wi 



Vi 



p(r) d^r / p(r') dV 



(9) 



We assume that Wy is proportional to the product of 



in appropriate units. And the total cost of running the 
network is proportional to the sum T -I- jZ with 7 > a 
constant that measures the relative importance of the two 
terms. Then the optimal network is the one minimizing 
this sum (25<. i26i| . 

Using again the conterminous United States as an ex- 
ample, we have first determined the optimal placement 
of p = 200 facilities which we then try to connect to- 
gether optimally. The number of edges in the network 
depends on the parameter 7. If 7 — > 0, the cost of 
travel Z vanishes and the optimal network is the one 
that simply minimizes the total length of edges. That 
is, it is the minimum spanning tree, with exactly p — 1 
edges between the p vertices. Conversely, if 7 — > 00 then 
Z dominates the optimization, regardless of the cost T 
of maintaining the network, so that the optimum is a 
fully connected network or clique with all ^p{p — 1) pos- 
sible edges present. For intermediate values of 7, finding 
the optimal network is a non-trivial combinatorial opti- 
mization problem, for which we can derive good, though 
usually not perfect, solutions using again the method of 
simulated annealing ^12]. 

There is, however, another complicating factor. In 
Eq. © we assumed that travel costs are proportional to 
geometric distances between vertices, which is a plausible 
starting point. In a road network, for example, the quick- 
est and cheapest route is usually not very different from 
the shortest route measured in kilometers. But in other 
networks travel costs can also depend on the number of 
legs in a journey. In an airline network, for instance, pas- 
sengers often spend a lot of time waiting for connecting 
fiights, so that they care both about the total distance 
they travel and the number of stopovers they have to 
make. Similarly, the total time required for an Internet 
packet to reach its destination depends on two factors. 



FIG. 5: Optimal networks for the population distribution of the United States with p = 200 vertices for different values of S 
and with 7 = 10~". 



the propagation delay proportional to the physical dis- 
tance between vertices (computers and routers) and the 
store and forward delays introduced by the routers, which 
grow with the number of intermediate vertices. 

To account for such situations, we generalize our def- 
inition of the length of an edge and assign to each edge 
an effective length 

hj^{l^6)k,+S (10) 

with < 6 < 1. The parameter 6 determines the user's 
preference for measuring distance in terms of kilometers 
or legs. Now we define the effective distance between 
two (not necessarily adjacent) vertices to be the sum of 
the effective lengths of all edges along a path between 
them, minimized over all paths. The travel cost is then 
proportional to the sum of all effective path lengths 

Z^J2^vh^ (11) 

i<j 

and the optimal network for given 7 and S is again the one 
that minimizes the total cost T + 7Z. (Since the second 
term in Eq. H1U|I is dimensionless, we normalize the length 
appearing in the first term by setting the average "crow 
flies" distance between a vertex and its nearest neighbor 
equal to one.) 

In Fig. I^lwe show the results of the application of this 
process to the lower 48 United States. When 6 = pas- 



sengers (or cargo shippers) care only about total kilome- 
ters traveled and the optimal network strongly resembles 
a network of roads, such as the US interstate network. 
As 6 increases the number of legs in a journey starts 
playing a more important role and the approximate sym- 
metry between the vertices is broken as the network be- 
gins to form hubs. Around 6 = 0.5 we see networks 
emerging that constitute a compromise between the con- 
venience of direct local connections and the efficiency of 
hubs, while hy d = 0.8 the network is dominated by a few 
large hubs in Philadelphia, Columbus, Chicago, Kansas 
City, and Atlanta that handle the bulk of the traffic. On 
the highly populated Californian coast, two smaller hubs 
around San Francisco and Los Angeles are visible. In 
the extreme case (5 = 1, where the user cares only about 
number of legs and not about distance at all, the network 
is dominated by a single central hub in Cincinnati, with 
a few smaller local hubs in other locations such as Los 
Angeles. 



V. CONCLUSIONS 

We have in this paper studied the problem of optimal 
facility location, also called the p-median problem, which 
consists of choosing positions for p facilities in geographic 
space such that the mean distance between a member 
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of the population and the nearest facihty is minimized. 
Analytic arguments indicate that the optimal density of 
facilities should be proportional to the population den- 
sity to the two-thirds power. We have confirmed this 
relation by solving the p-median problem numerically 
and projecting the facility locations on density-equalizing 
maps. We have also considered the design of optimal net- 
works to connect our facilities together. Given optimally 
located facilities, we have searched numerically for the 
network configuration that minimizes the sum of main- 
tenance and travel costs. A simple two-parameter model 
allows us to take different user preferences into account. 
The model gives us intuition about a number of situations 



of practical interest, such as the design of transportation 
networks, parcel delivery services, and the Internet back- 
bone. 
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